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Abstract 

Large earthquakes often trigger viscoelastic adjustment for years to decades depending on the 
rheological properties and the nature and spatial extent of coseismic stress. The 2006 My8.3 
thrust and 2007 M,8.1 normal fault earthquakes of the central Kuril Islands resulted in 
significant postseismic gravity change in GRACE but without a discernible coseismic gravity 
change. The gravity increase of ~4 uGal, observed consistently from various GRACE solutions 
around the epicentral area during 2007-2015, is interpreted as resulting from gradual seafloor 
uplift by ~6 cm produced by postseismic relaxation. The GRACE data are best fit with a model 
of 25-35 km for the elastic thickness and ~10'* Pa s for the Maxwell viscosity of the 
asthenosphere. The large measurable postseismic gravity change (greater than coseismic 
change) emphasizes the importance of viscoelastic relaxation in understanding tectonic 


deformation and fault-locking scenarios in the Kuril subduction zone. 


1. Introduction 

On 15 November 2006, the shallow part of the central Kuril Islands megathrust (depths <25 km) 
failed in a great My 8.3 earthquake following a swarm of thrust-faulting foreshocks in late 
September 2006 [Ammon et al., 2008; Steblov et al, 2008; Lay et al., 2009]. Following this event, 
distributed thrust-faulting aftershocks occurred along the forearc; additionally, a parallel band of 
primarily normal faulting events occurred in the outer rise region. Based on the predicted 
evolution of stress redistribution following the 2006 thrust event, Ogata and Toda [2010] 
hypothesized that continuous triggered slip occurred on the future fault plane of the January 2007 
event. On 13 January 2007, the great My 8.1 normal faulting earthquake ruptured a steeply 


dipping fault (depths of 4-33 km) beneath the trench slope [Lay et al., 2009]. Various studies of 
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seismic, geodetic, and tsunami data analysis report moment release between 1.6-5.3x10°' Nm 
for the 2006 event and 1.4-2.8x10*! N m for the 2007 event [Steblov et al., 2008; Lay et al., 
2009]. 

Earthquakes with moment magnitudes of My8.1—8.3 are nearly at, or below, the detection 
threshold by Gravity Recovery And Climate Experiment (GRACE). Furthermore, the coseismic 
gravity changes were dominantly a negative anomaly (due to crustal dilatation) following the 
2006 thrust rupture that was then disrupted by a positive anomaly (crustal compression) after the 
2007 normal faulting. We calculated a predicted total coseismic gravity change of <1 Gal at 
GRACE’s spatial resolution. This was certainly below the detection threshold [e.g., Han et al., 
2013]. However, following these events, the GRACE Release-05 Level-2 (L2) solutions 
[Bettadpur et al., 2014] exhibited an increasing gravity trend localized around the epicenters of 
the 2006 and 2007 earthquakes. This positive trend is consistently observed in various other 
GRACE solutions processed by different groups (Supporting information). It can be 
conveniently illustrated by using the CNES/GRGS web tools (http://www.thegraceplotter.com). 

Compared to the postseismic gravity change after the 2004 Sumatra-Andaman, 2011 
Tohoku-Oki and 2012 Wharton Basin earthquakes [e.g., Han et al., 2008, 2014, 2015], gravity 
change associated with the 2006-2007 Kuril events is distinctly different; namely, the magnitude 
of postseismic gravity change becomes by far larger than the coseismic change within a few 
years. We seek to explain such observations by testing alternate viscoelastic models of the 
asthenosphere beneath an elastic lithosphere. In this study, we address (7) how the 2006 M,8.3 
thrust and 2007 M,8.1 normal faulting earthquakes interplay to produce the observed gravity 
change, (i7) how the rheology of the asthenosphere plays a role in producing postseismic gravity 


changes, and (iii) how GRACE constrains the rheological structure in the central Kuril Islands. 
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We will present the details of GRACE observations and viscoelastic modeling results with a 
range of rheological parameters (thickness of lithosphere and asthenosphere viscosity) and of 
coseismic moment magnitude ratios, and compare with the previous Kuril GPS studies of Kogan 


et al. [2011, 2013]. 


2. GRACE observation of postseismic gravity change 

We used monthly GRACE gravity data from CSR L2 solutions presenting an effective resolution 
of 500 km (maximum spherical harmonic degree of 40). Although the GRACE RLOS5 L2 data 
are provided up to the maximum degree of 60, the resolution to ~300 km is applicable only at the 
polar region. The monthly gridded data, computed from the truncated spherical harmonic 
synthesis, were analyzed by applying the least squares fit of the parameters including annual and 
semi-annual sinusoids, bias, and a Heaviside step and linear trend after year 2007. Figure la 
presents the linear trend estimated between January 2007 and April 2015. The coseismic and 
postseismic effects after 2011 Tohoku-Oki earthquake dominate the trend map (Figure 1a). 
However, a statistically significant positive anomaly of 0.5 wGal/year is found around the 
epicenters of the Kuril earthquakes. 

We examined the finite fault models by Ji [2006, 2007] with the southeast dipping plane 
for the 2007 normal fault to compute the gravity effect as shown in Figures 1b-d. The synthetic 
calculation of coseismic (elastic) deformation after the 2006 earthquake results in the negative 
gravity change northwest of the trench due to crustal dilatation while a positive gravity (due to 
the hanging wall uplift) is compensated for by the ocean [e.g., Broerse et al., 2011] (Figure 1b). 
This is a typical pattern in coseismic gravity change after thrust earthquakes [Han et al., 2013]. 


Broerse et al. [2014] found the horizontal motion of the bathymetry is important to understand 
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the cosesimic gravity change. However, the horizontal postseismic motions during the 
observation period are ~2 orders of magnitude less than the coseismic horizontal motions, and 
hence they contribute a much smaller fraction of the net gravity change. The 2007 normal fault 
yielded instead positive gravity due to compression pushed by the downward motion of the 
hanging wall (Figure lc). Again, the gravity change due to density change is greater than the 
gravity change associated with the elastic vertical motion at this spatial scale. The gravimetric 
signals produced by the two earthquakes compensate each other and resulted in a gravity 
anomaly smaller than each of the gravimetric signals in the individual earthquakes (Figure 1d). 
The elastic gravity change is within +0.5 Gal and appears away from the epicenters at the 
spatial scale of 500 km. The solutions by Lay et al. [2009] including the southeast and northwest 
dipping planes for the 2007 rupture produced slightly different synthetics, but the above 
conclusion does not change (Figure S1). Elastic deformation alone does not replicate the peaked 
positive trend near the epicenter in the GRACE data. 

The time-series of the gravity change around the peak GRACE signal in Figure la 
(averaged within a spherical cap with a radius of 2.5° centered at 46.5°N, 154.5°E) is shown in 
Figure 2. The GRACE data (Figure 2, black dotted) presents scatter of 1.8 uGal (root mean 
square) with respect to the fit (green solid) as described above. The linear trend and annual 
sinusoid components are statistically significant while the step function and semi-annual 
components are not. The trend prior to the 2006 earthquake is not statistically different from 
zero (based on the t-test). The background gravity variation induced by other geophysical 
processes such as non-tidal atmosphere (6-hourly ECMWF operational data) and ocean mass 
(OMCT) [Dobslaw et al., 2013] is also shown (magenta solid), which was removed a priori from 


GRACE gravity solutions. The amplitude of seasonal atmospheric and oceanic variations is 
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estimated to be ~1 uGal. The GRACE data presents the un-/mis-modeled seasonal change with a 
amplitude of ~0.5 Gal (green solid). Ignoring the two-month time difference between the 2006 
and 2007 ruptures, the linear trend after year 2007 is estimated to be 0.45+0.08 pGal/year. The 
postseismic gravity trend estimate was consistent with the result based on the other GRACE 
solutions by CNES/GRGS, processed with a different inversion method and different atmosphere 
and ocean models (Figure S2) and by other groups (Figures S3 and S4). 

The gravity increase by ~4 uGal over 8 years from GRACE implies ~6 cm of the seafloor 
uplift around the area of ~500° km’ near the epicenters (0.7 uGal of gravity change per | cm of 
seafloor uplift - the Bouguer correction with the density contrast between the crust and ocean). 
Kogan et al. [2011 and 2013] presented a postseismic uplift rate up to 1.2 cm/year (KURILNET 
site MATC) from near-field GPS stations along the volcanic arc, which confirmed the 
postseismic deformation is caused predominantly by viscoelastic relaxation (rather than afterslip). 
We henceforth assume that the GRACE gravity changes are shaped primarily by viscoelastic 


relaxation and ignore the possible contribution of afterslip. 


3. Viscoelastic relaxation modeling 

For the Kuril earthquake sequence, the postseismic gravity change has emerged as the GRACE 
time series have lengthened even though the measured, combined coseismic gravity change was 
insignificant. In this Section, we examined in what circumstances the postseismic gravity 
changes could be larger than the cosesimic gravity change within a few years after the rupture by 
quantifying the gravity response functions (Green’s functions) of the viscoelastic Earth. We 
delineate the overall range of viscoelastic parameters important to control the postseismic gravity 


change by using a simplified point earthquake source and testing various synthetic structures. 
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Following Han et al. [2013], the gravity change by an earthquake point source is 


expressed as 





2( 0.0.1) =|) oe (@01)+( Moo) Bods) + (My) in!) 
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where a dislocation source is assumed to be at the north pole and described with the moment 


tensor components, M, M,,, M 0 Mo —M yg; and, M, je The respective Green’s function is 


Lure Suro Sure? Sueo-meg? MA Yiygs- They are dependent on the Earth’s viscoelastic structure 


and the depth of dislocation source. The first one g,,,. is constructed with the m=0 (monopole 


source) eigenfunctions of viscoelastic normal modes [Pollitz, 1997; Pollitz et al., 2006]. The 


pair of g,,. and g,,,, 1s determined by the m=1 (dipole source) eigenfunctions. The other pair 


Su0-moo AA Lig, 18 computed by the m=2 (quadrupole source) eigenfunctions. Therefore, there 


are in principle only three independent functions that govern temporal changes of gravity and 
displacement fields for each source depth with a given viscoelastic structure. 

Using the normal mode approach, we computed the viscoelastic (postseismic) and elastic 
(coseismic) Green’s functions. The seismic source was assumed to be at 20 km depth within the 
elastic layer. The viscosity for the upper mantle (220-670 km) and for the lower mantle (670- 
2891) was fixed to be 10°” Pa s and 10°! Pa s, respectively. The most sensitive parameters to 
govern temporal evolution of gravity are the asthenosphere viscosity and the lithosphere (elastic) 
thickness. Three different elastic thicknesses, 25, 40, and 115 km, were tested; below these 
depths the viscoelastic asthenosphere extends down to 220 km. Two different Maxwell 
viscosities, 10'* and 10'° Pa s, for the asthenosphere were tested (Table 1 of Supporting 


Information). 
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Figure 3 shows the ratio of the postseismic gravity response 5 years after the rupture 
relative to the coseismic change at each spherical harmonic degree (/) for the m=0, m=1, and 
m=2 responses, computed from six different cases for the viscoelastic structure. Note that such 
ratio is not sensitive to the seismic source parameters (except depth), but only to the Earth 
structure. Generally speaking, the shallower asthenosphere and the lower viscosity yield larger 
changes in the ratios. The ratio of the m=1 response is significantly smaller than the other 
response functions for all cases, which implies that the viscoelastic relaxation is governed more 
effectively by m=0 and m=2 patterns. We found that the asthenosphere viscosity of 10° Pas 
does not produce the postseismic gravity change that is significantly larger than the coseismic 
change regardless of the thickness of the overlying elastic layer. However, the asthenosphere 
viscosity of 10'* Pa s and an elastic layer thickness of <40 km can produce an order-of- 
magnitude larger postseismic gravity change than the coseismic gravity in 5 years. We extend 


the viscoelastic modeling to constrain such parameters in the following Section. 


4. Constraints on the asthenosphere structure 

We exploit the GRACE gravity observation to deduce the optimal parameters that may represent 
the rheological structure in the central Kuril Islands region. First, it is necessary to fix the 
seismic sources of the 2006 and 2007 earthquakes. Lay et al. [2009] indicates a large variability 
in seismic and geodetic source models (point sources as well as finite faults shown in their 
Tables 1 and 2). The moment release from the 2006 rupture has been estimated to be Mp=1.6— 
5.3x107' Nm including the Global Centroid Moment Tensor (GCMT) Myg=3.5x107! N m [Lay et 
al., 2009]. For the 2007 normal fault, it is Mj=1.4-2.8x107' N m with GCMT M,=1.8x107' N m. 


The moment release ratio of 2006 to 2007 ruptures varies from 1.1 to 3.8. Considering such 
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large uncertainty in the moment release from the two events, we decided to adjust the moment 
release ratio in addition to the viscoelastic structure parameters. Due to limited detection 
threshold of the GRACE measurements, it is not possible to constrain the moment release for 
each earthquake separately. 

The synthetic gravity change was computed using the finite fault models of Ji [2006, 
2007] and the software VISCOID for viscoelastic deformation [Pollitz, 1997]. Figure 4 shows 
the average rate of postseismic gravity change over 2007-2015 for the exemplary case of an 
elastic thickness of 32 km and the asthenosphere Maxwell viscosity of 10'* Pa s; the deeper (> 
220 km) upper mantle and lower mantle are as described in Section 3. The 2006 fault yields the 
positive gravity change primarily due to vertical motion of various density interfaces (such as 
seafloor and Moho) because viscoelastic deformation is shown to be nearly incompressible [e.g., 
Han et al., 2014, 2015]. The top of the asthenosphere is subject to the boundary condition of 
broad-scale shortening, which triggers postseismic uplift [Ve/osh, 1983]. In contrast, the 2007 
normal fault (for both northwest and southeast dipping planes) yields negative gravity change 
associated with postseismic subsidence under the top boundary of the asthenosphere subject to 
extension (See the schematic diagrams in Figure 4). The actual 2007 rupture, along with the pre- 
slip on this normal faulting rupture plane, in the outer rise caused contrasting vertical 
deformation to that triggered by the 2006 earthquake at the plate interface. 

We write the total postseismic gravity change, that is comparable to GRACE observation, 


from a given viscoelastic model as follows: 
g(t) = (DO) Sr (d,,V,t)+(1+ B)- 5007 (d,, v,t) ’ (2) 
where the synthetic gravity changes g,,,, and g,. , are evaluated at time ¢ with the viscosity V 


and the elastic thickness d, from the 2006 and 2007 ruptures, respectively. The unknown scale 
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factors, a@ and 2, are introduced to accommodate the uncertainty in the moment magnitude in 
each seismic model, respectively. They are not constrained independently but in terms of a ratio 


_(i+a) ms 
“(A 





where M a and M a is the scalar moment for the 2006 and 2007 


earthquake determined from a priori seismic models, e.g., Ji [2006, 2007], respectively. The 
ratio implies the relative size of the moment release of the 2006 to the 2007 rupture. We tested 
these three parameters within a reasonable range; i.e., d,=26-88 km, v=10'7-10" Pa s, and 7 
=1.0-4.0. 

We computed numerous cases of synthetic gravity changes and compared with the 
GRACE data. We quantified the goodness of the fit by computing the R’ statistics (i.e., variance 
reduction). Figures 5a-c illustrate the R’ statistics with variable viscosities and moment ratios, 
respectively, for d, = 26, 32, and 38 km. The results of d,>40 km were ruled out since the fits 
were significantly degraded within the specified range of V and y. The best fits are found 
consistently when the viscosity is close to 10'* Pa s. The larger elastic layer thickness d, 


requires the larger moment ratio Y (The identical results can be drawn based on the reduced chi- 


squared y’ statistics). The long-period seismic estimates of the relative moment ratio between 
the 2006 and 2007 ruptures are within 2—3 [7. Lay, 2015, personal communication]. With this 
seismic constraint, the case of d,=38 km can also be ruled out. The GRACE data are best fit 
with a model of 25-35 km for the elastic thickness and ~10'* Pa s for the Maxwell viscosity of 
the asthenosphere. 

The gravity time-series from three best-fit model parameters (color crosses in Figure Sa- 


c) were compared with the GRACE time series (Figure 5d). All cases reasonably fit the data 
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yielding R’ >0.75 and X ~1.0. We also tested the biviscous (Burgers) rheology that implements 
the transient (Kelvin) element in series with the steady-state (Maxwell) element [e.g., Han et al., 
2014, 2015]. The identical shear modulus was used for the Kelvin and Maxwell elements. We 
found that the transient viscosity is not as important as the steady-state viscosity to explain the 
GRACE data in 2007-2015. Any biviscous model with a Maxwell viscosity of 10'* Pa s and a 
Kelvin viscosity of 1-10x10'’ Pa s fits as well as the Maxwell rheology model with 10'* Pa s 
(Figure 5d). The Burgers or Maxwell model with the steady-state viscosity of 10'” Pa s, however, 
significantly underpredicts the gravity change (<1 uGal in 8 years) and does not agree with the 
GRACE measurement (~0.5 uGal/year). The spatial pattern of the average rate of the synthetic 


gravity changes in 2007-2015 is shown in Figure Se using the model parameters of d,=32 km, 
v=10'* Pa s, and 7=3.0. The primary positive gravity around the epicenters is consistent with 


the GRACE data in terms of the spatial location, the anomaly magnitude, and the temporal 


pattern, while the secondary smaller anomaly of the negative change is less discernable. 


5. Discussion 

Based on the analysis of 8 years of continuous gravity data, we estimated a Maxwell steady-state 
viscosity of 10'* Pa s for the asthenosphere, agreeing with Kogan et al. [2013]. Based on their 
analysis of 4 years postseismic GPS displacement data, Kogan et al. [2013] conjectured the 
Maxwell asthenospheric viscosity is lower than that in Chile and Alaska (10' Pa s), but these 
events were much larger in terms of spatial extent and moment magnitude. With the additional 
seismic constraints on the relative moment between two ruptures, our gravity analysis indicates 
the elastic layer thickness confined to be <35 km, which may be largely influenced by the mantle 


wedge beneath the active Kuril volcanic arc. Our preferred elastic thickness is smaller than 62 
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km fixed (not adjusted) by Kogan et al. [2011, 2013] for their GPS analysis, and rather 
consistent with the oceanic lithosphere thickness (30-40 km) from the flexural analysis of the 
bathymetry around the Kuril trench [Hanks, 1971] and with the numerical modeling beneath 
adjacent northeastern Japan [Muto et al., 2013 and Hu et al., 2014]. 

Lateral variations in rheology are likely around subduction zones, as differences in elastic 
plate thickness, elastic moduli and viscosity structure are expected between the landward and 
oceanic side of the subducting slab (e.g., Billen and Gurnis, [2003]; Abers et al., [2006]). This 
variation and the presence of the slab can substantially modify predicted postseismic deformation 
(e.g., Cohen et al., [1996]; Pollitz et al., [2008]; Hu and Wang, [2012]) and presumably the 
gravity signal as well. 

In our modeling, we neglected the lateral heterogeneity of the Earth’s structure that may 
bias our estimates of the elastic layer thickness. For example, the 2007 outer rise fault occurred 
in a thicker elastic layer [Lay et al., 2009] and is expected to have produced smaller viscoelastic 
relaxation than in a thinner elastic layer, which would lead to less disruption of the viscoelastic 
uplift produced by the 2006 rupture. However, the Maxwell asthenosphere viscosity estimate is 
consistent regardless of the assigned thickness of the overlying elastic lithosphere. As depicted 
in Figure 4, the 2006 rupture triggers an initially large viscoelastic uplift while the 2007 rupture 
adds the extensional boundary condition that pulls down the triggered uplift. What GRACE 
observed is the remnant of those contrasting postseismic vertical motions. 

As the gravity time-series lengthen, the sensitivity increases more drastically for 
detecting the postseismic trend than the coseismic change. Due to relocking of the shallow 
segment of the plate interface, a small part (as yet not discernible) of the postseismic gravity 


change will be constant interseismic deformation as well. While the magnitude of coseismic 


12 


266 


267 


268 


269 


270 


271 


272 


273 


274 


2Z1o 


276 


ZIT 


278 


279 


280 


281 


282 


283 


284 


285 


286 


287 


288 


gravity change is proportional to moment magnitude, postseismic (viscoelastic) gravity change is 
also dependent on the rheological structures that are variable from region to region. Our analysis 
of the My8.3 2006 and M,8.1 2007 Kuril earthquakes prompts reexamination of the GRACE 
data particularly with regard to postseismic gravity changes associated with even smaller 
earthquakes including the My8.1 2009 Samoa-Tonga earthquake doublet. This will become 
increasingly relevant as the GRACE time-series lengthens and the time series of gravity is 
extended with GRACE Follow-On starting in 2017 [Flechtner, et al., 2014]. 

This 2006-2007 Kuril sequence is one of five earthquake sequences that has produced 
persistent viscoelastic relaxation signals in GRACE gravity; the other cases are those following 
the 2004, 2010, 2011, and 2012 great earthquakes [e.g., Han et al., 2014, 2015, Broerse et al., 
2015]. The prevalence of viscoelastic relaxation was also observed in the seafloor data above the 
2011 Tohoku-Oki rupture [Sun et al., 2014]. The continuing analysis of the Earth’s rheological 
structure is imperative to understanding of the complete earthquake cycle processes. Wang et al. 
[2013] emphasized the importance of viscoelastic relaxation in drawing an adequate picture of 
fault-locking distributions at megathrust boundaries, noting that incorrect assessment of seismic 
and tsunami hazards may result from not considering the Earth’s rheological response. In 
addition, understanding of the Earth’s rheological response is also crucial for estimating climate 
signals (water and ice mass change) as the elastic and viscoelastic (un-)loading effects that exist 


in geodetic data such as GPS and GRACE need to be removed [J/vins et al., 2011]. 
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Figure 1. (a) The spatial pattern of postseismic gravity change from the linear trend of GRACE 
gravity data between 2007 and 2015 [wGal/year]. The sizeable 2011 Tohoku-Oki earthquake 
effect can be seen, southwest of the Kuril region. In this example, the GRACE solutions from 
CSR were used up to degree and order 40 equivalent to 500 km in spatial resolution; (b) and (c) 
The synthetic coseismic gravity change [Gal] after the 2006 thrust rupture and 2007 normal 
fault, respectively, computed using the seismic finite fault models and the isotropic PREM 
elastic structure; (d) The combined coseismic gravity change [Gal]. The calculated coseismic 


(elastic) gravity change alone does not account for the observed overall GRACE signals. 


Figure 2. The time series of the GRACE data around the peak GRACE signal in Figure la 
(within a spherical cap with a radius of 2.5° centered at 46.5°N, 154.5°E). The coseismic 
correction of the 2011 Tohoku-Oki earthquake (-0.5 uGal) was applied in the GRACE time- 
series following Han et al. [2014]. The GRACE data is shown in black with its error estimate, 
the data fit in green, and the a priori model variation of non-tidal atmosphere and ocean mass 
variability in magenta. The postseismic trends are by far larger than the background atmosphere 


and ocean mass variability. 


Figure 3. The ratio of the Green’s (eigen) functions between postseismic (viscoelastic) and 
coseismic (elastic) gravity change in 5 years after the rupture, respectively, in (a), (b), and (c), 
for the response functions of the monopole (m=0), dipole (m=1), and quadrupole (m=2) sources. 
The ratios are depicted for different spherical harmonic degrees up to 50. An earthquake point 
source is assumed at depth 20 km. For computation of viscoelastic normal mode relaxation, 


three different thicknesses of the elastic layer and two different viscosities of the underlying 
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asthenosphere were tested. In the monopole and quadrupole gravity responses, the postseismic 
gravity change was up to one order of magnitude greater than the coseismic gravity at the 


GRACE spatial scale, depending on the rheological structures. 


Figure 4. (a) The synthetic viscoelastic gravity change of the 2006 thrust fault in terms of the 
linear trend between 2007 and 2015 [wGal/year]; (b) The same but for the 2007 normal fault. 
The Maxwell asthenosphere viscosity of 10'* Pa s and the elastic thickness of 32 km were used 
in this example. Schematic illustrations of the positive (and negative) gravity rate due to 
postseismic uplift (and subsidence) under different boundary conditions of the thrust (and 


normal) faulting imposed on the top of the asthenosphere. 


Figure 5. (a)-(c) The R’ statistics between the GRACE gravity data and various synthetics. The 
R’ value of 0.7 means the variance of the GRACE data are reduced by 70% by removing the 
model time-series from the GRACE data. Three different model parameters were adjusted; 


elastic thickness d,, relative moment release Y, and asthenosphere viscosity V; (d) The time- 


series of the GRACE data, as in Figure 2, and the best-fit Maxwell rheology models (as marked 
with the color crosses in 5a-c) for the elastic thicknesses of 26, 32, and 38 km. Two examples of 
the Burgers rheology model with the steady-state viscosity of 10'* Pa s and the elastic thickness 
of 32 km and with the transient viscosity of 10'’ Pa s (Burgers!) and 10'° Pa s (Burgers2) are 
also shown. The synthetics with the thickness >40 km deviate significantly from the GRACE 
data and thus were ruled out and are not shown; (e) The spatial pattern of the synthetic 
viscoelastic gravity change from the best-fit model parameters, as marked with the green cross in 


(b). This is comparable to the GRACE data shown in Figures la. 
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